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NOTATIONS 


A a specified crack depth 

C constant in Paris's equation for crack growth 
H height of casing 
K stress intensity factor 

critical stress intensity factor 
stress intensity factor after N uses 
N number of uses of the motor case 
p Pressure 

R outer radius of the casing 
; o 

a surface crack depth 

critical crack depth 
a^ crack depth after N uses 

c q half the length of a surface crack 

c^ payload cost per pound 

c 2 cost of total payload 

c^ cost of articles and accessories at proof test 
c^ component cost 

component cost 
c^££ component cost 
c^ v component cost 
c,j, total cost 

f probability density function 

h constant in Paris's equation 

to initial case thickness 

p proof load factor 

t thickness of the case 

t„ thickness of the case after N uses 
x random variable representing crack depth 

z standard normal variable 

y,T) shape parameter 

At thickness decreased during grit blasting 
AK stress intensity range 

maximum initial crack depth, scale parameter 


iv 


V) " • 

a Q minimum crack depth 

(f) shape parameter 

y density of the material of the casing 

y yield stress 
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INTRODUCTION 


Structural components of a solid rocket motor case are considered to 
be! fracture critical whenever the game plan is to recover and reuse the 
moitor case for a designated number of missions. Proof tests, conducted on 
the case between missions, are also significant to rendering the 
structural components fracture critical. Proof load levels may 
significantly affect the design life of the structure. A fracture control 
plan is, therefore, necessary and is considered in the design of the 
case. 

In particular, this paper is concerned with the fracture control of 
the most critical membrane areas of the case. All discussions and 
methodologies presented in this paper can, however, be used whenever 
similar fracture critical structures of a reusable space vehicle system 
are designed. Some modification might be necessary in particular 
structures. Significant loads are applied to the motor esse during flight 
and water recovery operation of each mission. The applied stresses from 
all other events during the mission are assumed in this analysis to be not 
significant enough to result in cyclic or time dependent crack growth. If 
the test or analysis indicate the possibility of other critical loading 
events they can be included in the fracture control plan by extending the 
reported analysis. Before each mission, the case is also subjected a 
proof test. The loads ' applied during the proof tests can result in 
significant amount of crack growth. Grit blasting is assumed to be used 
between each mission. This reduces the effective depth of cracks and the 
thickness of the membrane by a selected amount. While the effective depth 
of crack is reduced, the refurbishment grit blasting operation has the 
effect of increasing the applied stresses. This necessitates a larger 
initial thickness of the membranes than that would be required otherwise 
Therefore, any design of the membrane of the case must arrive at an 
initial wall thickness t, the thickness t that will be decreased between 
each mission and the proof load factor p. For example, a large value of 
initial wall thickness results in increased reliability, but results in 



the need for increased propellant, increased cost of operation and reduced 
pay load capability. On the other hand, a small initial wall thickness 
increases the probability of failure and the resulting loss of the reusable 
space vehicle system and the pay load. Therefore, there is a need for 
optimizing the initial wall thickness. Similar arguments can be presented 
to explain the need for selecting the other design variables such as At and 
p by optimizing the desired objective function of cost and weight. 

In general, these design variables depend on the probability 
distribution for the initial flaw sizes present in the membrane, applied 
stresses during the use of the vehicle, crack growth characteristics of the 
material, fracture control plans, specified reliability bounds, weight and 
cost considerations. The paper describes a reliability-based procedure 
that can be used to select the design variables of a solid rocket motor 
case in a reusable space vehicle system by using probabilistic fracture 
mechanics and cost or weight considerations. 

Method of Approach 

It is assumed that careful nondestructive inspection (NDI) 
techniques can detect initial cracks greater than the surface length of 2c q 
and depth of a Q with 100% success. Sometimes, it is assumed that 

cracks corresponding to surface length 2c = 0.1 inch can be identified 100% 
of the time.l If the corresponding maximum depth is 0.05 inch there is no 

possibility of existence any initial cracks of depth larger than 0.05 inch. 

! 

Such an initial crack depth distribution is assumed to be analytically 

j 2 . . 

represented by Johnson S^ distribution. Reasons for this assumption can 

be explained as follows. One of the requirements of any assumed 

distribution is that the minimum and maximum crack depths be bounded within 

finite limits. Depending on the thickness and the available techniques of 

non destructive inspection techniques, there is a finite maximum depth of 

possible crack. It is not infinity as is provided by distributions such as 

normal distribution, gamma or log-normal distributions. The minimum value 

of depth of crack can be assumed to be zero or a small number . Such a 

distribution can be obtained as the transformation of the usual normal 

variate. One such transformation is the following. 
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e ^ x £ e + X (1) 


In this equation, z is the standard normal variable and x is the variable 
of interest i.e., the crack depth. The four available parameters are y , 
T),i e, & X . The minimum and maximum crack depths fix e & X respectively. 
The parameters can be called shape parameters and can be determined from 
percentiles of the observed data. 


The density function for the probabilistic model is written as 
follows 


f / a \ = -11 £ 

ao' o' 2 tt (a *■ e) (X - a + e) 
o o 


exp \-k 
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This empirical distribution is called Johnson distribution. It should 
be noted that it is possible to obtain other empirical distributions to 
represent the crack depths. 
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This probability distribution for initial crack depth changes after 

each mission, each proof test and each time the material as removed from 

the wall thickness. The change in distribution after each mission and 

each proof test is due to the crack growth resulting from the applied 

stresses. This crack growth also depends on the present length of the 

crack, applied stress and the material that are responsible for the crack 

growth, in this analysis, the applied stresses and material properties are 

assumed to be known deterministically. If the initial crack length were 

also known deterministically the crack length after each use can be 

3 .4 

determined from equations such as Paris' equation , Foreman's equation 
or Collipriests equations'*. Because initial crack lengths are not known 
deterministically, crack length after each use of the vehicle is again 
another probabilistic distribution that has to be estimated. 


is 
uses 


The cumulative density function (CDF) for crack length after N uses 


denoted by F(a^). 


This represents the probability that a^^A after N 


Each use is defined as one flight, one proof test and a material 
removal. Crack growth due to time related effects such as stress 
corrosion have been neglected. 


If F(a^) is known, the probability distribution for the stress 
intensity factor can be obtained from the knowledge of the applied 
stresses. The probability distribution F(K^) for stress intensity factor 
can be used to estimate the probability failure which is the 
probability of stress intensity factor K N greater than or. equal to the 
critical stress intensity factor during the projected design life of the 
structure. The critical stress intensity factor is denoted by K £ , In 
this analysis, stresses and the material properties are assumed to be 
known deterministically. However, the applied stress changes after each 
use due to material removal. Therefore, the probability of failure can be 
expressed as the probability of a ^a^. In this expression a c is the 
critical crack depth that can be obtained from the critical stress 
intensity factor and the applied stress corresponding to that particular 
mission. This relationship between the stress intensity and the applied 
stress is discussed in the next section. 



Stress Intensity Factor 


For the analysis of the stress intensity factor in the membrane, an 
infinite plate model with elliptical surface flaws that are oriented per- 
pendicular to the applied stress has been assumed. The relationship 
between the stress intensity factor, the applied tensile stress and crack 
depth is given by^ 


(3) 


where 

2 

Q (|) = 0 2 - 0.212 (4) 


In this equation, 0y is the yield stress and 0 is a function of the 

2 

ratio of crack depth to crack length (a/c). Variation 0 with (a/c) is 
given in reference 1. 

Because the crack depth a is a random variable the stress intensity 
factor K is also a random variable. In general, both crack depth a and 
crack length 2c are random variables and there is a need for a joint 
distribution for a and c. In this analysis, only the crack depth is 
considered as the random variable. It is also assumed that the 
probability distribution for crack depth a is known initially and is given 
by a Johnson distribution . The density function for the distribution 
: :f --- 5 




is given in equation (1). This probability distribution for crack depth 
changes with use. The next step will be to determine the change and the 
new probability distribution after each flight and proof test. 

Probability Distributions for Crack Depth After Use 

The following symbols are used to properly account for the changes in 
probability distributions. 

f(a Q ): Probability density function for the 

initial crack depth 

F(a ): Cumulative distribution function for 

o 

initial crack depth 

F(a ): Cumulative distribution function for 

°P 

initial crack depth after the first proof test 
F(a„): Cumulative distribution function after 

N flights and (N+l) tests 

F(a„ ): Cumulative distribution function after 

Np 

N flights and N proof tests. 

F(a„): Cumulative distribution function after 

N 

material removal from the wall thickness. 

3 

Then, using Paris's crack growth equation, 

i = c W n (5) 


where C and n are empirical constants. Alternately, the rate of crack 
growth can be assumed to be given by Foreman's equation^ or Collipriest ' s 
equation"*, if they are found to represent the situation more accurately. 
For example, Collipriest ' s equation can be written as follows: 


da 

dN 


™ ■ ex P 




In K - lnAK 0 „ ,-l llnAK- % (InK (1-R) + InAK )l 

c u tanh / c p ' , 


\ (In K (1-R) - InAKo 


+ in |c exp ln V lnAt! n ) 1] 


( 6 ) 
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where Paris's constants are used. By integrating either of the selected 
equations (5) or (6) crack depth after N+l uses can be determined if the 
crack depth after N uses and N proof tests are known deterministically, 
i.e. , 



Similarly, crack depth after the proof test can be determined from 
equation (5) or(6) if the crack depth before proof test is known 

deterministically, i.e., 


a NP " a NP | a N| 


( 8 ) 


These functions represented by equations (7) or (8) can be determined 
analytically or in the form of quadratures from equation (5) or (6). From 
equation (7), a N+ ^ can be obtained for every known value of a Np . 
Similarly, a^ p can be obtained for every known value of a N from equation 
(8). However, both a„„ and a„ are random variables in t e present 
analysis. In this case, equation (7) can be used to obtain the 
probability distribution for a^ + ^ if the probability distribution for a^ p 
is known by using the principle of transformation of random variables. It 
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should be noted that all equations similar to (7) or (8) involving crack 
depths are increasing functions. This property is useful in transforming 
the random variables. 

For example, the probability density function for a N+ ^ can be written 
as follows 



(9a) 


similarly 


£ <V ' 


[> M 


da 


N 


da. 


Npj 


(9b) 


Equations (7) and (8) can be written for every value of N from zero 
to the projected number of uses. 

Details of obtaining these equations for the membrane of the solid 
rocket motor case, with the expression for stress intensity given by 
equation (2) and Paris' equations for crack growth, is discussed in the 
Appendix I. 

The next step is to obtain a tool for change of probability 
distribution due to the material removal from the wall thickness. 

Material Removal and the 
Change of Probability Distribution 

Due to material removal after each use, the effective crack depth is 
reduced by A t . Thus, new crack de^th is 


8 


^ ' 


At 


( 10 ) 


It is assumed that At is a constant. Thus, by using the principles of 

2 

transformation of random variables , the probability density function for 

a„ can be written as follows. 

N 

P(\) = f + At) (11) 

In this equation, p ( a n ) represents the density function for a.^ and f 
represents the functional form of the probability density function for a^. 

Probability of Failure 

By following the method discussed in the preced.ing two sections 
probability density function for crack depth can be obtained after every 
flight, proof test and material removal. From the density function, 
cumulative probabilities can be obtained by integration. Integration 
after the transformation of variables as discussed in equations (9), (10), 
and (11) needs the determination of appropriate limits of integration 
consistent with the transformation of variables. This is also discussed 
in the Appendix I. If F(a^) represents the cumulative density function 
after N flights & N proof tests the probability of failure is given by the 
probability of a > a c g- The quanity of a^ corresponds to K £ and the 
applied stress at the use. 
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It is to be noted that the probability of failure changes with 
different selections of the initial wall thickness t, increased loading 
due to proof test, the material removed At and the number of designated 
number of missions. The increased loading due to proof tests is denoted 
by a factor p. A cost function or a weight function can be formulated from 
this knowledge of probability of failure and other related unit-cost or 
weight. Such a cost or weight function depends on t, p, and number of 
missions N. It is possible to select these design variables by minimizing 
the cost or weight function subject to appropriate realiability bounds. 
The effect of nondestructive inspection (NDI) is indirectly related to 
initial flaw distribution. Additional NDI effects such as the rejection 
of structures are not considered in the analysis. However, they can be 
included as units related to the probability of failure. A numerical 
example is illustrated in the next section to illustrate the developments 
of the paper. 

Numerical Example 

For the numerical example, it is assumed that the Johnson S K 
distribution for the initial crack depth is such that the minimum crack 
depth is zero and the maximum crack depth is 0.1 inch. Paris's .^equation 
for crack growth is assumed with 


C' = 0.847 x 10" 

n = 3.0 


2 

The variation of 0 with (a/c) as shown in figure 1 is approximated by a 
quadratic relation. 

The primary objective of reusing the solid rocket motor case is to 
reduce the cost of operation of the reusable space vehicle system in which 
it, is used. However, as the number of uses is increased, the probability 
of failure increases because of the propagation of the crack depth. On 
the otherhand, smaller number of uses increases reliability and also the 
cost is distributed over a smaller number of uses. This means the casing 
has to be replaced after a fewer number of uses. 
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A larger initial thickness would increase the weight of the casing 
and costs more in terms of payload. But the probability failure is less 
if the thickness is more. The proof test factor p and the material 
erosion At are kept constant in this example. However, they also can be 
varied and their effect on total cost can be considered in the most 
general case. The total cost function c^,, therefore, comprises the 
following component costs. 

i) Initial cost of the casing, c^, 

ii) Expected cost of flight failure c^, 

iii) Expected cost of proof test failure and 

iv) cost due to multiple usage, c . 

The initial cost c^ is given by the product of the weight of the casing 
and the cost per pound of the system, i.e., 

c i - " < 2 V„ - S 2 ) » V (12) 

where R ^ outer radius of the casing 

t = thickness of the casing at the 

„th . 

N cycle 

H= height of the casing 
y= density of the material 
c = payload cost per pound 

The expected cost of flight failure is the product of the probability of 
flight failure and the entire payload cost, i.e. 

C ii “ P N ' C 2 (13) 
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fch 

where is the probability of failure at the N flighty is the total 
cost of the payload. Similarly the cost of proof test failure is 


c . . . = P„ c 
ill Np 3 


( 14 ) 


th 

where p is the probability of failure at the NT proof test and is 
cost of articles and accessories of proof test. Finally, the cost due to 
multiple usage is given as follows: 


0.3 

(N) (15) 


Thus, substituting all the components, the total cost function c^, is given 
by the following equation 



c_ = c. + c. . + c... + c, 

T i ii in iv 


(16) 


The following numerical values are used^’^ in evaluating equation (16). 


y = 0.3 lbs/cubic inch 

H = 816 inches 

R = 72.5 inches 

o 

= $1624 per lbs. 

C 2 = $250 x 10 6 

C 3 = $2xl0 6 

Results 


The Initial thickness t D is varied from 0.535 inch to 0.435 inch in steps 
of 0.005 inch. Also, 1 % of the initial thickness is eroded after each flight 
The total cost function is calculated for various initial thicknesses and 
use cycles by means of a digital computer. Figure 2 illustrates the 
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variation of the cost function with t and N. If is obvious that as the 

o 

number of uses increases, the minimum occurs at a higher initial 
thickness. For example, for 18 missions the minimum cost occurs at an 
initial thickness of 0.48 inch. The initial thickness to give minimum 
cost for 20 mission cycles increases to 0.497 inch, for 22 missions the 
thickness required is 0.512 inch. 

Figure 3 delineates the variation of reliability with initial 
thickness, after 20 missions cycles. The reliability corresponding to the 
minimum cost for 20 uses is 99.3%. If this reliability is not ade- 
quate, then a higher initial thickness should be used even though the 
total cost will be higher than the minimum. 


General Procedure 

Based on the preceding example, a general procedure can be delineated 
in the following steps. 

2 

1, Obtain the parameters of the Johnson distribution for 
the initial flaw size. 

2. Obtain the stress in the membrane from the known geometry of 
the case and wall thickness. 



In the equation p is the proof stress factor. During flight, 

p is replaced by a value of 1. Pressure P is the MEOP pressure 

on the case and R is the radius of the case, 
s o 

3. Obtain the new CDF and density function for the crack depth after 
the proof test. 
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4. Obtain the new CDF for the crack depth during the flight 
following the proof test. 

5. Estimate the probability of failure. 

6. Compute the cost function parameters. 

7. Obtain the new CDF after the material removal. 

8. Repeat steps 2 to 7 for the new thickness and the next 
mission until the total number of missions are complete. 

9. Change t and N and repeat the calculations as necessary. 

10. Select the design variables for the minimum value of the 
objective function subject to reliability constraints. 

A computer program has been written to carry out these steps .{see Appendix 

I* 

Conclusions and Recommendations 

This paper has demonstrated that the reliability analysis based on 
probabilistic fracture mechanics can be used to optimize the selection of 
the design variables of a solid rocket motor case. In particular, basic 
design variables such as the thickness and projected design life as well 
as: the fracture control variables such as the proof factor and material 
erosion can be included in the analysis. Accuracy in estimation of the 
initial flaw size distribution is reflected in the assessment of the risks 
involved in the design. By knowing the risks involved in the design, 
weight and cost can be reduced from those obtained by the conventional 
deterministic analysis and use of arbitrary safety margins. 

This report is only a first step in the development of procedures 
based probabilistic fracture mechanics. Additional work that is 
necessary can be listed as follows: 
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1. A more accurate analysis can be obtained by considering the 
joint distribution for the crack depth and crack length 

along the surface. 

2. Accurate methods of estimation of the probability distribution 
for the initial flaw size distribution should be developed. 

3. In particular, effects of water impact and time dependent 
crack growth, stress corrosion, should be considered. This 
is particularly important if the missions are spaced over 
years . 

4. Uncertainties in external loads and material properties 
should be considered. 

5. Accuracy of the different models for crack growth (in the 
point of view of probabilistic fracture mechanics) should 
be evaluated. 

6. Alternate fracture control plans and more accurate stress 
intensity measures based on cylindrical geometry can be 
considered. 

7. Cost of NDI efforts in relation to the cost that will be 
incurred by additional safety factor should be evaluated 
in the point of view of improved reliability. 


8. Thermal effects should be considered. 
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APPENDIX I 


This appendix describes estimation of the new CDF of crack depth after use 
from a knowledge of the old CDF and probability density before use. 

Crack Growth Rate 

The rate at which the crack depth increases is given by Paris's 
equation as follows. 


~ - C(AK) n . 0.847;rl(f l8 (AK) n 


For subsequent convenience in algebra, the value of n is taken to be 3.0. 
By substituting for c and n, 


da 

dN 




a 

%r 

= 0.847 

e *i 

m. 

,a N 2 

[ c 5 + c 2f + \ 

m 


X10 


-18 


(Al) 


Simplifying this further, 


da 

dN 


1.5 


* c. 


c 5 + c 2 


<!> 


,a ‘ 
c 3 ( c> 


(A2) 


where 
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c, = 0.847 x c.^ x 10 ^ 
6 4 


(A3) 


Separating the variables a and N in it follows that 


dN 



i C 5 + C 2 ( c ) + C 3 ( c ) 


1.5 

da 


(A4) 


Integrating both sides between state (1) and state (2) the following 
equation is obtained 



(A5) 


In order to evaluate the integral on the right hand side, it is found 
necessary to expand the numerator of the integrant binomially. 

Now consider the numerator of the integrand with C, = 1. Neglecting 

3 5 

terms of higher order than (a/ c ) « it follows that 


18 



Letting 


and 


[i + c 2 (i) + c 3 (f) 2 ; ! 


1.0 + 1.5 c. (r) 
2 c 


+ 1.5c 3 + 1.5 (0.25)] 


❖ 


+ [o.75c 2 c 3 - 0.25(0.5) 2 c 2 3 ] (|) 


(A6) 


P 1 = c 1,5 C 2 


(A7) 


*2 [ 

r. V 


2 2 ) 1,5 C 3 + 0,375 c 2 

c 


(A8) 


P 3 “ ~3 i0.75c 2 c 3 - (0.25) 2 c 2 3 | 


(A9) 


Then, it follows that 


I 1 


+ c 2 (f) + c 3 (f) ! = 1.0 + P 3 a + P 2 a 2 + P^ 3 (AW) 


Substituting in the integral the following result is obtained 
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1 

0.5 


, v-0.5 , 

(a) + 


0.5 


•(a) 


0.5 



(a) 


1.5 



Solution of a ^ as a function of 


Substituting the limits of integration in A(ll) 


c 6 <H 2 - »j) - -u 2 -°- 5 + 2?^ 2 ) 0 - 5 + |y 2 


? s "0.5 n c ic ? 5 

+ | P 3 a 2 + 2a i ' 2P l a l " | P 2 a 2 ‘ | P 3 a i 

Rearranging and neglecting terms of order higher than three, : 
the following equation 


(a ) + p (a ) + q(a ) + 7 = 0 


where 


1.0 

P ‘ | P 1 P 2 • 4 P 3 


(4P. 


8P 9 ) 
3 1 
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(All) 


.5 

(A12) 

.t reduces to 
(A13) 

(A14) 



- 1.0 

q ■ 8 P.P, - 8 P 
3 L L 5 


( 8P 1 + c x 2 ) 


and 


y • 


8 p 1 p 2 

3 1 * 



Now, the three roots of this cubic equation, (a ^) 1 are given as 


(1) 


= A + B - P 


1 3 


(2) A + B A - B 

1 ~ Z + z 

\£T 

(3) A + B A - B 

1 z z 

V^T 


where 

A.J/-b + Jb 2 + iJ 
J 2 y U 27 

B - J- i ^ + i 

yj 2 \ U 27 

a = O (3q - p 2 ) , b » 1_ (2p 3 - qpq + 27y) 

27 


(A15) 

(A16) 

follows 

(A17) 

(A18) 
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Transformation 


Probability density of a ^ is given by 


% < V ■ ^ % < V 


(A19) 


CDF of a 2 is then 


a ( a o) 

1 


r a 2 

c 

I f (a ) da 

= 1 f 

/ - a 2 2 

2 / a 

J al 2 

* o 

W 


f £ a, < a l> da l ' 

\ <a i>’ 

Jo 1 

1 


( a 1 ) d a L 


VV 


(A20) 


(A21) 


2 

where F (a,) is the CDF of Johnson S_ distribution . 
a^ i d 

Now, it is needed to obtain a^ as a function of a , No. of cycles etc. 
This can be done by solving the polynomial equation obtained previously in 
terms of a^ and treating a^, N^ and as constants. The infinite degree 
polynomial equation is truncated at the 3rd degree for convenience. 

Of the three roots only one will be the real root because of the 
physical nature of the problem, r :y a^Ca^) 


Then by substituting in the expression for the CDF of a 0 




I (a2> jo f&1 


(a L ) da x 


(A22) 
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Substituting the values 

for point (i), 


Substituting the values 

for point (ii) 



1.0 + c 2 (0.5) + c^ (0.25) = 1.5 

(A27) 

or 

2 c 2 c 3 = 2.0 

(A28) 

Substituting the values 

for point (iii) 


i 

1.0 + c 2 + c^ = 2.5 

(429) 

or 

^2 + ^3 =1 * 5 

(430) 

Solving equations (2) and (3) simultaneously 



c 2 = 0.5 

(431) 

and 

c 1.0 

(432) 


Thus the chosen parabolic fit is as follows 


0 2 = 1.0 + 0.5~ + ~2 (433) 

c 
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By hypothes is , the initial flaw a, has a Johnson distribution 
Ai so . there is a functional relationship between the initial flaw size s. 


and the suoseouent rlaw size 


"ter a. cycles. 




a random variable because a^ is a randan variable by hypothesis. 
Having known the range space of a T the range space of e 9 can be derived 
from the functional relationship between a ? and a... 'Thus, if the lower 
limit of a ^ is zero, it follows from the functional relationship between 
a^ and that the lower limit of a ? is also zero. Next, if Che upper 
limit of a^ is a^,, the corresponding upper limit for a 9 can be obtained by- 


solving the cubic relation between a., and a„ , as a function of the number 
2 11 
of cycles N = “ N^. 
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i M MAIN 


74/74 CPT=1 


FTN - • E + - ? 3 


?■■/ 12/ 0- . 


PROGRAM MAIN ( I NPUT, OUTPUT , TA Pr3 = lNPUT-rTAP- PsOUTPUT ) — 

**•*************■■*••¥•♦**. * *♦♦** **•*•*■¥• ***♦*■**♦.*♦♦ * W****¥*«4f4 *■,**.**.* + ..*, 

THIS PROGRAM CALCULA T£S Th£ RELIABILITY OF THE’ S-M CASING AT TH _ 
EACH PROOF TEST AND 'JS£ CYCLE GIVEN TH:. INITIAL F l. A N OIST-I5UT I CN 
VARIOUS INITIAL T HIC<NtSS^E AND LAUNCH NUMR-RS. 

DIMENSION F(30.1> .Y(.5G,1> ,ZZ1 ( r ) 

COMPLEX A-A.BE.A1HAT 

READ (5,*) ALAMOA.EATA . GAMA 

£PS=0.G 

WRITE (6*100) ALAMOA* -ATA, GAMA *E°S 
ENCFILE6 

FORMAT <1H1,/,5X, "PARAMETERS OF THE JOHNSON SB DISTRIBUTION ARE 


C-AMA =“,FE.2 


PS XL ON =' 


= ** ♦ FI 0 . .5 ) 

ONE PROOF 


END OF EACH CYCLE. 


*-**■» **•.*-**• 


1ILAMBDA =*' *F«*.1.2X.* , -ATA =" . F5 . 2 • 2-X « "C-AMA =“ *F~ .2 . 2X-, M PS ILON = ", 
1F5.2) 

DC 5000 1=1,21 
AKP=1.05 

•*AKP" IS THE PROOF-TEST FACTOR - 

WRITE(6,102)AKF 

FORMAT (1H1,/,5X. M PROOF - TEST FACTOR =*\F15.3> 
STHIK=0.486*(1.-(I-11)/100.) 

--WRITE (6,-10 3) STHIK - 

‘FORMAT ( / , 5X » *' INITIAL THICKNESS OF SRM CASE =",F10.3) 

DO 4000 N=l, 25 

M N" IS THc. NUMBER OF THE CYCLE CONSISTING CF ONE PROOF-T^ST AND ON,. U: 

DO 3 C 1-0 L = 1 » 2 - - — - — - 

IF(L.RQ.l) WRITE (6,104) N 
IF(L.tQ.2> WRITE (6,10 5) N, 

FORMAT (///,5X, "PROOF TEST NUMBER =",Il3) 

- FORMAT (5X, "LAUNCH NUMBER OF- SRM-GASE -",17) — 

THIS SECTION CALCULATES THE THICKNESS AT THE END OF EACH CYCLE, 

THIK = STHIK- (( (STHIK/ 100 , >*(N-1) ) > 

WRITE (6, 106) THIK - - - 

♦2S*lIi*^*n5J***»*i»2**.f*2**tf **"i*i*;*t*« *.♦♦♦**. *«*•*«»******«**»*»* 

THIS SECTION CALCULATES THE APPLIED STRESS (SIGM°) 

- SIGMP=AKP*950.7*72.5/THIK ... 

ACR= ( <93 500/ SIGMP) **2.0 ) / (1.2*3.143) 

IF (L.EQ.2) S1GMP=SIGMP/AKP 
WRITE (6,108) SIGMP 

**♦**♦*♦***♦ **,»***.*«,;****** ***.*..*, 

THIS SECTION CONSIDcRS THE CUBIC APPROXIMATION. 

C = 0.1- 

Cl=1.0 - * 

C2=0.5 
C3=l . 0 

C4=(SQRT (1.2*3. 1-*?)) ‘SIGMP 
■C6 = 0.847MC4**3,G) 

C6=1.0E-16*C6 

P1=1.5*C2/C 

P2*f (1.5*C3)* (1.5*0. 25*C2*C2) )/ ( C*C ) 

»3=(<1.5*C.5*C2*C3)-<0.0625*(C2**3.0>>) /<C**3.04 

A2=ACR 

Zl = 2 , C *P3* (A2**2 . 5)/6 .0 
Z2 =2 .0 *P2* (A 2**1 . 5)/ 6 ,0 
Z3=2.0*P1*(SQRT(A2)) 


O O 

*11 

•55 

g£? 


i 


112 


co 

o 


C 

C 


llH 


C 

C 


c!5.v) 


Z4 = * (2.0/SQST (42) )-(Co*N) 

C7*Z1+Z2+Z3+Z4 
Zl*(L.0*Pl*Pll-(8.0*P2/3. 0) 
Z2 s (6.0*Pl*P2/3.0)-(8.0*P3/5.Q) 

P=Z1/Z2 

Zl*-((6.0* o l)+(C7*C7) ) 

Q=Z1/Z2 
R=4.0/Z2 

ABAR=( (3.0*Q) - (P+P)) /3.0 
Z1=2.0MP**3.0) 

Z2--9 *P*Q 
Z3=27*P 

BMZ1 + Z2+Z3) /27.C 
Zl=(B*B/4.0) + ( (ABA R/ 3.0)* *3.0) 

FORMAT (/," 0 1 S C"R E M I N _ N I SQUARE 
AA= ((-B/2.0)+CSQRT(CMPLX(Z 1,0.0) > ) 

8B= ( ( *B/2 • 0 ) -CSQR T(CMPLX(Z 1,0.0) ) ) 

Ql=REAL(AA) 

Q2=AIMAG (AA) 

RR=(SQRT( (Q1**2.0)+(Q2**2.Q) ) )** ( 1 . 0/3 . 0) 
THET=ATAN(Q2/Q1) 

AR1=RR*COS(TH£T/3.0) 

AIl*RR*SIN(THET/3.0) 

AR2=RR.*C0S( (THET+ (2.*3.l47) 1/3.0) 

- AI2=RR*3IN ( (THf_T-M2. *3.14.7)) /3.0) 
AR3=RR*COS((THET+(4.*3.147))/3.0) 

AI3=RR*SIN( (THET+(4.*3.147> >/3.0) 

Ql=REAl(BB> 

-02=AIMAG (88 ) 

RR=(SQFT ((Q1**2.0)+(Q2**2.0) ) )** (1. 0/3.0) 

THE T=ATAN(Q2/Ql) 

9R1=RR*C0S ( THE T/3 • 0) 

• - BI1=R.R*SIN (THET/3.0) • - *- 

BR2=RR*C0S( (THET+ ( 2. * 3 . 14 7 ) ) /3 . 0 ) 

BI2=RR*SIN( (THET+ (2. *3 .147) ) /3.0) 

BR3=RR*C0S ( (THET+ (4. *3 . 1*- 7 ) ) /3 . 0 > 
-8I3=PR*SIN((THET+(4.*3,1^7))/3,0) 

A1HAT=AR1+8R1- (P/3.0 ) 

AA=REAL (A1HAT) 

*«******¥********' 

THIS SECTION CONSIDERS THE QUADRATIC APPROXIMATION. 

Z1=(4.0*P1*P1) -(8.0* “2/3. 0) 

Z2=-((8.0*P1)+(C7*C7) ) 

P=Z 2/Z1 
Q-4.0/Z1 

A 1 = ( -P+SQRT ( (P*P) -(4 . 0*Q> ) ) /2 .0 
WRITE (6, 114) A1 

FORMAT (5X, “UFFE R LIMIT OF &i = *',.2>.o) 

A1HA1=A1 

IF (A1HA1.GE.ALAMOA) Y(N1.1)=1.0 
IF (A1HA1 .GE . ALAMOA) GO TO 300: 

T H 13 SECTION C AL CUlA Tcb THE RELIABILITY 

N1=1E 

DO 3000 K=2.N1 
AJ.=A1HA1* (K-l) /F l JAT ( M i-1) 

Z1=E AT A* ALA MO A /SORT (2 .0*3. 1-7) 
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Z2*1»0/(A1* (ALAMDA-A1) » 

Z4=EATA*AL0G(A1/(ALAMCA-A1> > 

Z3=EXF(-((GAMA+Z4>**2.0)/2.Q) 

FCK,1»=Z1*Z2*Z3 

3000 CONTINUE - 

F ( 1* 1 ) =0 • 0 

CALL INTGRl(1,A1HA1,F.Y,N1> 

3005 WRIT£(6tll6>Y(M,l) 

116 -FORMAT 15 X f M P*06 ABItl T V OF NO FAIlUR. - ' - - - - 

ZZl(L)sl.-YCNl.l) 

3010 CONTINUE 

C ***#♦•♦•***♦.#***********♦*****♦*♦.**#**♦#*♦******♦**#*♦*■+**■*.**■**♦*♦*♦****** 

•C THIS- SECTION- CALCULATES THE TOTAL COST - FUNCTION . - - „ 

cii=(i .o /float (nm»o .30 
Cl=162c. 

C2-250.E+6 

G3=2.0E + 6 - - — - 

WT*3. 143*616. * 0 . 3M C 145 . * THIK> - f THIK*TH IK )) . 

CTOT=(WT*Cl» +(ZZl<2>*C2)+(ZZi (1 > *C3 > + <C3* Cl 1 > 

WRITE (6,119)CT0T 

118 TORN A T ( 5 X * "-* * ***** *** *********** TOTAL-COST FUNC T IGN-= F-l 5 . 5T 

4000 continue 
5000 CONTINUE 
999 STOP 



OF POOR qual: 


INt INTGRL 


’4/74 0PT = 1 


FTN -t« 6 + -2 i 


SUBROUTINE INTGRL. (M,T,f,-*SN) 

01 HE NS ION A(3Q,3Q),B(30,30),C<3Q,30) ,F<30,1) ,Y(3Q,1) 
H = T/ (720 .* (N-l) > 

00 20 1 = 1 , N 

DO -20 J=1.N — 

B ( I , J ) =0 ■ 0 
A(I.J)=0.0 
20 . CONTINUE 

DO 24 K = 1,N - 

OO 24 J= 1 ,K 
B(K,J>= 1.0 

24 CONTINUE 

A ( 2 , 1 ) =251. - - — ' 

A ( 2 , 2 ) = 646 . 

A(2,3)=-264. 

A < 2 » 4 ) =106. 

A (2, 5) =-19. - 

A (N,N-4> =-19. - ' 

A (N,N-3) =106. 

A (N»N-2> =-264. 

A (N,N-1) =646. ~ • ■ - 

A (N,N) =251 . 

A (N-l , N-4) = 11 . 

A (N-l.N-3) =-74. 

A ( N-l , N-2) =456. 

A( N-l, N-l) =346 • 

A ( N-l ,.N) =-19 . 

J = N-2 

OO 25 1 = 3, J 

A (I, 1-2) =-19 .' 

A (I, 1-1) =346. , 

A (I, I) =456. 

A ( 1,1 + 1) =-74. - - - 

A.CI, 1 + 2) = 11 , 

25 CONTINUE 

DO 30 1=1, N 

•OO 30 J=l~, H— — — - - 

C ( I, J ) =0 . 0 
DO 30 K=1 * N 

C(I.J) =G(I, J)+(B(I,K) )*(A(K,J)*H) 

,30 CONTINUE 

OO 35 1=1, N 

y (i, i)=o.o 

DO 40 J=1 » N 

40 Y (I,1)=Y (I,1) + (4C*(I,^))*(F( J,1>) ) 

35 . CONTINUE 

IF (M.EQ.l) GO TO 45 

M = M -1 

OO 65 K=1,M - - - r- - 

„ DO 50 L= 1 , N 

50 F(L,1) =Y (L.l) 

DO 55 1=1, N 

Y (I, 1 ) =0 *0 - - . - 

DO 60 J=1*N 

60 Y(I,1)=Y(I,1)+(C(I,J)*F(J,D) 

55 CONTINUE 

65 CONTINUE - 


INE INTGRL 


74/74 OPT =1 


F r N - . c + - 2 


45 - DO 80 1=1, N • 

80 Y ( I • 1 ) = A BS ( Y ( 1 , 1 1 ) 

RETURN 
ENO 
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Vc 


Figure No. 1. Variation of Shape Factor with a/c 




33 



• t> 

inickness (in) 


ij'jufo 2 . Mini! vi Sa t i on '"urves 




